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CHAPTER I 


INTRODUCTION AND LITERATURE REVIEW 


I.A Overview 

Accurate life assessment of structural components may require advanced life 
prediction criteria and methodologies. Structural components often exhibit several 
different types of defects, among the most prevalent being surface cracks. 

A semi-elliptical surface crack subjected to monotonic loading will exhibit stable 
crack growth until the crack has reached a critical size, at which the crack loses stability 
and fracture ensues (Newman, 2000). The shape and geometry of the flaw are among the 
most influential factors. When considering simpler crack configurations, such as a 
through-the-thickness crack, a three-dimensional (3D) geometry may be modeled under 
the approximation of two-dimensional (2D) plane stress or plane strain. The more 
complex surface crack is typically modeled numerically with the Finite Element Method 
(FEM). A semi-elliptical surface crack is illustrated in Figure 1-1. 
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Figure 1-1 Surface Crack in a Plate 

Characterizing surface crack growth and fracture under monotonic loading 
requires knowledge of the material behavior and stress state surrounding the crack front. 
In cases where the plastic zone surrounding the crack tip is of small magnitude relative to 
the distance to the nearest boundary, Linear Elastic Fracture Mechanics (LEFM) may be 
applied for a simple determination of failure loads and related quantities. High levels of 
plasticity may necessitate use of Elastic-Plastic Fracture Mechanics (EPFM). 

A mathematical description of the conditions required to induce crack growth in 
an elastic body was first presented by Griffith (Griffith, 1920) in the form of an energy 
balance equation. Crack-tip stress field expansions for an elastic body were later derived 
by Irwin (Irwin, 1956) and Williams (Williams, 1957). The crack-tip stress field 
expansion is dominated by a constant within the first term, the stress intensity factor 
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(SIF) K as denoted by Irwin. Similarly, the path independent ./-integral was proposed by 
Rice (Rice, 1968) as the dominant parameter in the elastic-plastic stress field expansion. 
Dodds et al (Dodds, 1993) remark that two fundamental concepts underlie both LEFM 
and EPFM: 

1) the relevant crack-tip singularity dominates over microstructurally significant 
size scales 

2) the parameter K or / uniquely scales the amplitude of the near tip field. 

I.B Finite Element Analysis 
I.B.l Overview 

The embedded elliptical crack in an elastic body was first studied by Irwin (Irwin, 
1962). Irwin provided the foundation for semi-elliptical surface crack Finite Element 
Analysis (FEA) investigations conducted by Ayres (Ayres, 1970) and Levy et al. (Levy, 
1971), who used 3D elastic-plastic small strain formulations to obtain the plastic-zone 
shape and stress distribution around the crack front. McMeeking and Parks (McMeeking, 
1979) and Shih and German (Shih, 198 1) later utilized FEA to examine components 
under applied monotonic tension and bending in order to characterize the evolving stress 
fields near the crack-tip. A useful summary of the advances in the characterization of 
elastic-plastic crack-tip fields is presented by Parks (Parks, 1992) 

Raju and Newman (Raju, 1979) (Newman, 1981) performed 3D elastic analyses 
to obtain K for semi-elliptical surface cracks for a range of crack sizes and loading types. 
The Raju-Newman K solutions have since been expanded by Fawaz and Andersson 
(Fawaz, 2004) who analyzed the corner crack at a hole configuration. The solutions were 
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extended to larger ranges of crack-depth-to-thickness and crack-depth-to-width ratios by 
utilizing higher levels of mesh refinement and large-scale computational resources not 
available to Raju and Newman. Trantina et al. (Trantina, 1983) also performed elastic- 
plastic surface crack FEA to establish the limitations of LEFM and to compute the J- 
integral for small cracks. Parks and Wang (Parks, 1992) presented /-integral values for 
surface cracks determined using detailed finite element solutions, and studied the effects 
of local crack front constraint on the fracture process. 

I.B.2 Finite Element Mesh Design 

One of the difficulties in applying FEA to surface cracked geometries lies in the 
generation of meshes suitable for accurate calculations in the crack-front region. The 
near tip stress fields in a cracked body are dominated by stress and strain gradients 
normal to the crack front due to the immobility of the material in front of the advancing 
crack. For reliable finite element calculations, the geometry must be adequately 
discretized in the region local to the crack front to capture these gradients. Faleskog 
(Royal Institute of Technology, Stockholm, Sweden) developed a code to generate semi- 
elliptical surface crack meshes based on a right-hand curvilinear elliptic coordinate 
system derived by Timoshenko (Timoshenko, 1970). Historically, meshes generated 
with the Faleskog code have produced reliable results (Faleskog, 1995) (Gao, 1998) 
(Aveline, 1999). However, it is difficult to generate a mesh without high aspect ratio 
elements. More recently, Structural Reliability Technologies of Boulder, CO (www.srt- 
boulder.com ) have developed a commercial software package FEA-Crack capable of 


5 


generating meshes for various cracked geometries, including surface cracks. The FEA- 
Crack surface crack mesh is based on a rectangular coordinate system, and utilizes a 
highly discretized tube of elements around the crack front. The mesh is generated using a 
proprietary code; a license must be purchased in order to use the software. A typical 
FEA-Crack surface mesh is presented in Figure 1-2. Meshes generated with FEA-Crack 
were used for all analyses reported herein. 



Figure l-2.a Typical Semi-Elliptical Surface Crack Finite Element Mesh 



Figure l-2.b Finite Element Mesh Details around the Crack Front 
I.C Fracture Prediction Methodologies 

The development of fracture prediction methodologies has been the topic of much 
research. The most accurate approach to predict structural integrity is to forecast the 
fracture process including initiation of crack growth, the extent of stable growth, and 
failure. Typically, fracture analyses are conducted by calculating a well-defined fracture 
mechanics parameter and comparing it to a critical value that has been detennined 
through material testing. Hult and McClintock (Hult, 1956) discovered that under 
conditions of limited crack-tip plastic deformation, the details of the local elastic-plastic 
fields could be uniquely related to a single macroscopic parameter, such as K or J, scaling 
the intensity of crack-tip deformation. Parks (Parks, 1992) comments, ”... that local 
crack tip fields can be characterized by a single parameter, and further, that fracture 
processes are driven by these fields, there exists a mechanistic rationale for constructing 
‘single parameter’ fracture mechanics correlations of crack extensions.” 
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I.C.l Single Parameter Fracture Criteria 

The amount of local yielding around the crack front dictates the appropriate 
parameter choice. In cases where the relative plastic zone size is small, as in brittle 
materials, K is the fracture controlling parameter. However, a larger plastic zone size 
indicates behavior typical of ductile materials, in which case the ./-integral has been 
proven as a suitable parameter. 

The SIF has been correlated to the energy release rate in a cracked body. It has 
been widely used to predict crack extension by characterizing failure as the point where K 
is equal to the plane strain fracture toughness, Kjc. Reuter et al (Reuter, 2002), compared 
the plane strain fracture toughness K IC and K p k (the peak K value around the perimeter of 
a part-through crack) for different materials loaded in tension and bending. Using the 
conventional criterion for mono tonic loading to failure, ratios of K p k/ Kic were found to 
be greater than 1.0 and in some cases were greater than 2.0, implying that conventional 
practices were conservative. 

The ./-integral has been used to correlate the initiation of crack growth in 
plastically deforming solids. When high levels of plastic defonnation are present, the 
relationship between the /-integral and the crack-tip stress field lose a direct correlation 
(McMeeking, 1979) (Shih, 1981). The loss of /-dominance signifies a loss of constraint 
in the body and lends support for the incorporation of a second parameter in the fracture 


criterion. 
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I.C.2 Two Parameter Fracture Criteria 

In an effort to more accurately predict failure under monotonic loading conditions, 

incorporation of a constraint term as a second parameter in fracture prediction criteria has 

been proposed by several investigators, including Hancock et al. (Hancock, 1991) 

(Hancock, 1993). Constraint refers to the buildup of stresses around a crack front due to 

restraint against in-plane and out-of-plane deformation. Newman et al. (Newman, 1995) 

present a precise description of constraint: 

Strain gradients that develop around a crack front cause the deformation in the 
local region to be constrained by the surrounding material. This constraint 
produces multi-axial stress states that complicate stress analyses and influence 
fatigue crack growth and fracture behavior. The level of constraint depends upon 
the crack configuration and crack location relative to external boundaries, the 
material thickness, the type and magnitude of loading, and the material stress 
strain properties. 

Constraint has often been used within fracture mechanics in a qualitative manner, 
such as plane-stress or plane-strain constraint. However, efforts to quantify the influence 
of constraint on fracture have been the subject of much recent work. In order to use 
constraint in fracture prediction, the crack front stress state must be resolved with a 
numerical parameter defining the level of constraint along the crack front. McMeeking 
and Parks (McMeeking, 1979) expressed the plastic stress concentration factor as a 
measure of constraint 


Kop = maxJ^rr-l (1-1) 

where <% is the normal (crack opening) stress, R ( ./_ is the lower yield point, and max x is 
the maximum quantity in the x-direction. Several researchers (Rice, 1969) (Hancock, 
1976) (McClintock, 1979) used a more general definition of constraint given by <j m / o vm , 
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where a m and <j vm are the mean and equivalent stresses in the neighborhood of the crack- 
tip, respectively. Sommer and Aurich (Sommer, 1991) analyzed the mean-stress-to- 
equivalent-stress ratio for surface cracked specimens and showed how constraint (as 
defined in this manner) affected stable crack-growth behavior under monotonic loading 
conditions. Hancock et al. (Hancock, 1991) proposed the use of T stress, the stress that is 
in-plane and parallel to the crack surfaces, while others such as O’Dowd and Shih 
(O’ Dowd, 1991) used the Q stress as a measure of stress triaxiality around the crack 
front. Newman et al. (Newman, 1993) proposed using the nonnal stress in the near-tip 
stress field as a measure of constraint. This measure of constraint was referred to as the 
global constraint factor a g , the average normal stress acting over the plastic region 
through the thickness of a through crack. The hyper-local constraint factor a/,, developed 
by Aveline and Daniewicz (Aveline 1999), is defined as the average of the normal-stress- 
to-flow-stress ratio along a line emanating from the crack front to the plastic zone 
boundary along the crack plane. Newman et al. (Newman, 1999) and Reuter et al. 
(Reuter, 2002) used ah to predict initiation and fracture of surface cracks in brittle 
materials under tension and bending loads and correlated these results with cracked 
through-the-thickness bend specimens. Fracture initiation was predicted to occur at the 
load corresponding to the maximum ah K value along the surface crack front. The 
approach predicted initiation load within ±20%, displaying the viability of a/, as a 
fracture criteria constraint parameter. Aveline and Daniewicz (Aveline, 1999) developed 
ah for a range of crack sizes and loading types for brittle materials, but the approach has 
not yet been applied to ductile materials. 
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The objective herein is to extend the above work for application to a ductile 
material. ,7- integral and constraint factor distributions for a range of surface crack sizes 
and loading conditions were obtained. The fracture initiation location is then predicted to 
occur at the point of highest 7a/,. Finite element analyses were conducted on a wide 
range of surface crack configurations, and the /-integral and a/, values were calculated 
for each model. To verify the validity of the developed fracture prediction model, the 
Idaho National Engineering and Environmental Laboratories (INEEL) tested a large 
number of surface crack specimens, each loaded monotonically to failure. These 
specimens and particular surface crack configurations were modeled with the FEM. The 
analytical and experimental data is presented and the potential of Ja/, as a fracture 
criterion is discussed. 



CHAPTER II 


NOMENCLATURE AND CONSTRAINT DEFINITION 


II.ACrack Geometry and Loading 

Cracked bodies are characterized by geometrical parameters describing the size 
and shape of the flaw. In the case of the surface crack, as shown in Figure 1-1, the 
notation is as follows and is presented in Figure 2- 1 : crack depth a, crack half-length c, 
specimen half-width w, specimen half-height h, and specimen thickness t. When 
discussing surface cracks, it is common to describe the size of the crack in terms of 
geometrical ratios relating the crack to the containing body, where the crack-depth-to- 
specimen-thickness ratio is a/t and the crack aspect ratio is a/c. Parameters that describe 
behavior along the surface-crack front require an angular measure defining the location of 
a point on the crack front; the preferred nomenclature is the parametric angle (/). 
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Figure 2-1 Surface Crack Configuration 


Considering a point on the crack front, (x- h yj, the parametric angle may be calculated 
using Eq. 2.1. The free surface of the crack is defined as the location where </>= 0, and the 
deepest point of penetration where (j) = n/2. 


(f> = sin 



( 2 - 1 ) 


\a J 

The surface cracks modeled in this investigation have been subjected either to an 
applied monotonic tensile load S corresponding to a tensile stress <Jt or to a uniform 
bending moment M with a corresponding maximum bending stress of magnitude <J B . 
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Surface cracked plates exhibit two planes of symmetry, one along the length ( y-z 
plane), and the other along the width (x-v plane). The planes of symmetry and coordinate 
system orientation are presented in Figure 2-2. 



Figure 2-2 General Finite Element Model and Applied Loading 

The presence of geometrical symmetry greatly reduces the modeling requirements and 
eases the computational burden. Also evident in Figure 2-2 is the crack plane, the x-v 
plane. The crack plane is the location of most concern for the analyses presented. The 
shaded region on the crack plane will be referred to as the uncracked ligament or the 
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material ahead of the crack front, while the white region on the crack plane will be 
denoted as the cracked ligament or the material behind the crack front. 

II.B Constraint Definition 

The definition of constraint considered is an extension of the global constraint 
factor a g presented by Newman et al. (Newman, 1993). When a cracked body is 
subjected to an applied load, a small region of material directly ahead of the crack front 
will be elevated to a stress level beyond the material yield strength <j 0 due to the 
immobility of the material ahead of the crack. Considering the von-Mises yield criterion, 
the material is said to have yielded when the equivalent stress or von-Mises stress <j vm has 
reached cr 0 . The von Mises stress is determined by the overall 3D stress state present in 
the body, see Eq. 2-2. 

Consider the simple 2D plane stress case of a notched body subjected to a tensile 
load, Figure 2-3. At an infinitesimal point A at the notch tip, the stress in the x-di recti on 
<j x must be zero to satisfy the free-surface boundary condition. For the point A to be 
yielded, the stress in the v-direction <j y must be equal to cr 0 . Moving a small distance 
ahead of the notch tip to point B, the free surface boundary condition is no longer present, 
thus <j x is no longer forced to a value of zero, and <j v may be elevated above <j 0 by a factor 
a. This factor defines the level of local normal stress constraint present in the plastic 



zone of the material. 
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Figure 2-3 Two-Dimensional Constraint Definition 

In three-dimensional stress space, the local constraint is magnified by the 
presence of stress triaxiality and the resistance to in-plane and out-of-plane deformation. 
The higher the constraint level the more plastic deformation near the notch will be 


restrained. 


CHAPTER III 


LABORATORY TESTING OF SURFACE CRACKS 


III. A Test Plan 

INEEL was contracted to perform mechanical testing on plates containing surface 
cracks by NASA and the FAA. A plan was presented to INEEL outlining the material 
selection, surface crack fabrication specifications, and testing guidelines. The material 
selected for testing was D6AC steel, heat-treated to a ductile condition. The heat 
treatment as dictated by the ATK-Thiokol standard was: 

1. Austenitize under controlled atmosphere at 1615°F+/-25°F for 2.5 hours 
minimum. 

2. Quench in molten salt bath at 325°F maximum initial temperature, 15 minutes 
minimum. 

3. Air cool to 175 °F 

4. Snap temper in molten salt for 3 hours minimum at 3 10°F to 345°F 

5. Clean metal to remove all salt 

6. Temper to meet mechanical property requirements. Minimum of two 
temper cycles shall be used. Tempering temperature is 1070°F - 1 1 15°F for 6 
to 7 hours. Cool components in air to 175°F max between tempering cycles. 

After heat treatment, the material must satisfy the following property requirements: 

1. Ultimate Tensile Stress: 200 ksi min - 225 ksi max 

2. Yield Stress: 180 ksi min 

3. Percent Elongation: 8% min 

4. Reduction in Area: 25% min 

5. Kic: 90 ksiVin 

The D6AC steel used for specimen fabrication satisfied the above requirements. 
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Twelve different surface crack shapes and sizes were identified as the crack 
configurations to be tested. The same configurations were to be tested monotonically 
under either remote tension or bending loads. The crack configurations are presented in 
Table 3-1: 


Table 3-1 Proposed Surface Crack Configurations 


BENDING LOADS 

TENSILE LOADS 

a/c 

a/t 

a/c 

a/t 

0.90 

0.20 

0.50 

0.75 

0.90 

0.20 

050 

0.75 

0.72 

0.20 

0.50 

0.75 

0.72 

0.20 

0.50 

0.75 

0.48 

0.20 

* 

0.75 

0.48 

* 

0.50 

0.75 

0.10 

0.20 

0.50 

0.75 

0.10 

0.20 

* 

0.75 


The crack sizes denoted by * show a negligible change in Ki around the perimeter of the 
surface crack as calculated by the Raju-Newman equations, thus were not proposed for 
inclusion in the test plan. The specimens were fabricated and a triangular starter notch 
was electrical-discharge machined (EDM) in the center of the specimen to aid crack 
initiation. All specimens, regardless of load type were pre-cracked under remote cyclic 
bending loads of unknown magnitude 1 to reach the desired initial crack size specification. 
Obtaining a specific crack size is a difficult procedure, thus the initial surface crack sizes 
were not identical to the proposed configurations. After pre-cracking, 14 tension and 8 
bending specimens were available for testing. The specimens were loaded monotonically 
until a 5% potential drop was recorded, indicating that a small amount of crack extension 
had occurred. The load at the 5% potential drop was recorded. After the first occurrence 

1 Pre-cracking for the initial crack formation was not performed at INEEL and the load levels were not 
recorded. Shear lip formation was not observed on the fracture surfaces, indicating the pre-cracking levels 
were likely within reason. 
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of crack extension, the specimens were cyclically loaded at a reduced load level to mark 
the location and extent of stable crack growth. The specimens were subjected to three 
instances of crack extension and cyclic marking before loading to complete failure; 
however, only the first instance is considered in this study. The specimen and crack 
dimensions as well as maximum applied stress at the onset of first crack extension are 
presented in Table 3-2. a (tension) and Table 3-2.b (bending) below. 



19 


Table 3-2.aTension Specimen Specifications 


Specimen ID 

a (mm) 

c (mm) 

t (mm) 

w (mm) 

<y T (MPa) 

AT-01 

1.22 

1.6 

6.26 

25.28 

1300.8 

AT-02 

1.46 

1.76 

6.32 

25.46 

1294.2 

AT-04 

3.20 

4.44 

6.35 

25.40 

1152.3 

AT-05 

3.33 

4.75 

6.35 

25.40 

1111.7 

AT-06 

3.27 

4.56 

6.34 

25.50 

1140.1 

AT-07 

4.13 

6.78 

6.40 

25.39 

1006.3 

AT-08 

4.27 

6.97 

6.33 

25.36 

992.5 

AT-09 

4.16 

6.63 

6.24 

25.49 

987.8 

BT-01 

4.61 

9.80 

6.35 

25.40 

841.0 

BT-04 

5.89 

23.31 

6.35 

25.40 

271.4 

CT-01 

3.49 

6.44 

6.22 

25.44 

872.5 

CT-02 

3.51 

6.65 

6.20 

25.37 

998.4 

CT-03 

3.38 

6.56 

6.29 

25.31 

1044.2 

DT-02 

1.73 

6.42 

6.35 

25.40 

1205.6 


Table 3-2.b Bending Specimen Specifications 


Specimen ID 

a (mm) 

c (mm) 

t (mm) 

w (mm) 

Ob (MPa) 

AB-01 

4.356 

7.076 

6.37 

25.375 

1720.6 

AB-04 

3.185 

4.543 

6.38 

25.375 

1574.5 

AB-07 

1.516 

1.822 

6.37 

25.350 

1969.6 

BB-01 

4.620 

9.860 

6.35 

25.385 

1369.8 

BB-04 

5.495 

23.460 

6.32 

25.375 

448.3 

CB-01 

1.712 

2.690 

6.36 

25.37 

2010.3 

DB-03 

2.314 

6.487 

6.36 

25.35 

1772.2 

DB-04 

3.880 

15.850 

6.34 

25.385 

1047.7 


After loading each specimen to failure, a high-resolution digital image was taken of the 
fracture surfaces to show surface crack pre-cracking shape and size and the amount of 
crack extension. An image showing the details of the fracture surface is given in Figure 
3-1. Sample tension and bending fracture surfaces are provided in Figures 3-2. a and 3- 
2.b, respectively. 
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Figure 3-1 Surface Crack Fracture Surface Details 



Figure 3-2.a Typical Surface Crack Fracture Surface (Tension) 
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Figure 3-2.b Typical Surface Crack Fracture Surface (Bending) 


In addition to the surface crack specimens tested, single-edge-bend (SEB) specimens 
were fabricated and tested to provide plots of /versus crack extension Aa. While not 
considered as part of this research, they may prove beneficial for future research. 




CHAPTER IV 


FINITE ELEMENT ANALYSIS OVERVIEW 

WARP3D release 1 5 (WARP3D, 2004) was used in analyzing the 22 surface 
crack models (14 tension, 8 bending). FEA-Crack version 2.5.625 was employed to 
generate the surface crack meshes and the WARP3D input files. A preliminary 
verification for /-integral calculations was performed to ensure that the WARP3D 
solution parameters were being used correctly and that the surface crack meshes were 
adequately refined. The material model, mesh characteristics, and solution parameters 
are presented. 

IV.A Preliminary Verification 

/-integral calculations were verified against those published by Parks (Parks, 
1992). Parks performed 3D elastic -plastic FEA on surface cracked plates under varying 
tensile and bending loads and calculated the /-integral as a function of (j> around the crack 
front for each load level. To verify the calculation of /using WARP3D, finite element 
models were constructed of identical crack size, material model, and applied load level to 
those of Parks. The / solutions obtained from WARP3D were plotted against the Parks 
solutions for both the tension and bending cases. The / values were normalized by 
(£ 0 <j 0 tZ~), where e 0 and <j 0 are the yield strain and yield stress, respectively, t is the 
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thickness, and 27 is a loading parameter given by the applied stress divided by the yield 
stress. The comparisons are shown in Figures 4- 1 .a and 4- 1 .b. While the J solutions 
compare within reason, differences remain evident. The lack of agreement is likely due 
to a combination of mesh refinement limitations within the Parks solution and different 
methods of calculating the ./-integral (Parks used the Virtual Crack Extension method 
whereas WARP3D uses the Domain Integral method). In view of these differences in 
analyses, WARP3D was considered a reliable means of calculating the /-integral. 



Figure 4-1. a /-integral Verification for Models Subjected to Tension 
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Figure 4- l.b /-integral Verification for Models Subjected to Bending 
IV.B Material 

Tensile tests were performed on the D6AC steel at INEEL on specimens cut in 
both the transverse and longitudinal rolling directions. The resulting material properties 
were essentially isotropic. The results from all tests were averaged to obtain a single 
value for use in the finite element analyses and are summarized in Table 4-1. 

Table 4-1 Thiokol D6AC Steel Material Properties 


Yield Stress, cr 0 

1329.7 MPa 

Ultimate Tensile Stress, a„ 

1434.0 MPa 

Young’s Modulus, E 

209.7 GPa 
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An incremental plasticity, Ramberg-Osgood (Ramberg, 1943) material model was 
utilized in the finite element simulations. The Ramberg-Osgood model is defined by Eq. 
4.1, where cr 0 is the reference stress (yield stress), e 0 is the corresponding reference strain, 
n is the hardening exponent, and k is the fitting constant. 



(4-1) 


Using this stress-strain relation, a Ramberg-Osgood stress strain curve was lit to the 
tensile test data. The input parameters used were <j 0 = 1329.7 MPa, s 0 = <j 0 / E, n = 50, 
and k = 0.3 15. The curve-fit is presented in Figure 4-3 with the average tensile test data. 
Excellent agreement was obtained. 
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Figure 4-2 Ramberg-Osgood Curve Fit to Tensile Test Data 
IV.C Analysis Specifications 

The Bdisop element type was used for the analyses. It is an 8-noded isoparametric 
element and employs a tri-linear displacement field (WARP3D, 2004). The WARP3D 
sparse solver was used in the analyses. Within the finite element code, the von-Mises 
yield criterion and its associated flow rule were used. Linear-kinematic element 
formulations (small strain) were used. 
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IV.D Boundary Conditions 

The finite element model was constrained in a manner to simulate the symmetry 
planes of the surface crack (as outlined in Chapter 2). The y-z plane of the model was 
constrained in the x-direction and the un-cracked ligament on the crack plane was 
constrained in the 2 -direction. To prevent model translation, a single node located at the 
farthest point from the model origin ((x, y, z) = (w, t, 0)) was constrained in the y- 
direction. 

IV.E Loading Specifications 

The finite element models were loaded on the far face at z = h, in the z-direction. 
For the tensile cases, a uniform stress ay was applied on the element faces at this 
location, and for the bending cases, a linearly varying stress with <j B at the upper surface 
(y = t) and -a B at the lower surface (y = 0) was applied. The applied stresses for each 
specimen were equivalent to those presented in Tables 3-2a and 3-2b. The stresses were 
applied incrementally in 40 equal load steps, from zero to the maximum applied stress, in 
order to aid solution convergence. 

IV.F Convergence Problems 

Of the 14 tension and 8 bending models selected for analyses, two tension (AT-02, BT- 
04) and one bending model (BB-04) would not converge on a solution due to high levels 
of plasticity. These models are not considered in the presentation of results. 



CHAPTER V 


RESULTS AND DISCUSSION 


V. A Constraint Calculations 

The hyper-local constraint factor was calculated as a function of the parametric angle 
around the crack front for each surface crack model. A mathematical expression for the 
constraint calculation is shown in Equation 5.1. 


a 0 = 




S(<f) 

f 

f \ 


J 

0 

UJ 


ds 


(5-1) 


Figure 5-1 provides a graphical representation of the path definition. 



A FORTRAN routine ALPHAH was developed to calculate a constraint distribution for 
each surface crack model. The source code for ALPHAH is presented in Appendix A. 

The stress for each node on the crack plane 
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was taken as an average of the surrounding Gauss point stresses, and used as input for the 
routine. The nodes possessing a von-Mises stress greater than or equal to cr 0 were used to 
define the plastic zone around the crack front. Each node lying on the crack front 
represents an individual Rvalue, as defined in Eq. 2.1, for which a constraint value is 
calculated. The yielded nodes which lie in a path perpendicular to the crack front S 
emanating from each crack front node are sorted into subsets, and average normal stress 
to flow stress ratio for each subset is considered as the a/, value for the corresponding (f). 
The hyper-local constraint factor distribution along the crack front is then plotted versus 
(j). 

V.B /-Integral Calculations 

FEA-Crack greatly simplified the calculation of the /-integral through automatic 
generation of WARP3D input files containing the appropriate / calculation commands. 
The domain integral method is used by WARP3D to calculate the /-integral (WARP3D, 
2004). The average /-integral values for each (f> location are output to a results file that is 
used to generate plots of /versus (f) for each model, /-integral calculations were 
nonnalized by the product of the applied stress ( u T or <j B ) and the thickness t. 

V.C Fracture Initiation Location 

The high-resolution fracture surface images provided by INEEL were analyzed to 
determine the location along the crack front that exhibited the largest amount of crack 
extension. This is an arduous task and is highly subject to the interpretation of the 
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analyst. Each fracture surface image was digitized, and the (f) location of maximum crack 
extension on both sides of the initial surface crack was visually selected and recorded. 
Figure 5-2“ shows a typical fracture surface and corresponding ^ locations. The two (j) 
values were then averaged to obtain a single location of maximum crack extension for 
each test specimen. The point of maximum crack extension is considered the fracture 
initiation location, and will be used to validate the developed fracture criterion. 



Figure 5-2 Typical Maximum Crack Extension Location Measurement 


1 Note that in the figure shown, the crack did not grow symmetrically around the EDM notch, implying that 
it may not have been in the center of the specimen. 
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V.D Discussion of Results 

Surface cracks under tension and bending possess different characteristic crack 
front stress fields, ./-integral and constraint calculations for tension and bending thus 
differ in characteristic shape. The shape of each curve is related to the plastic zone size 
around the crack front. The ./-integral signifies the extent of stress and strain elevation 
along the crack front, so it shares a direct relationship with the amount of crack extension 
present. Conversely, increased plasticity signifies a loss of constraint; hence, constraint 
and plastic zone size exhibit an indirect relationship. The calculated a/, and ./-integral 
distributions were plotted for each model as well as the product /«/, normalized in the 
same fashion as the /-integral. A comprehensive set of plots for each crack configuration 
and loading type is presented in Appendix B. 

V.D.l Tension Specimen /-integral and ai, Distributions 

The /-integral and constraint variation around a monotonically tensile loaded surface 
crack front share the same typical shape. The plastic zone of a surface crack under 
tension appears as a bulge just beneath the free surface, and decreases to a constant value 
as the deepest point of penetration is approached. The /-integral and a/, distributions 
mimic this as both display a steep gradient just below the free surface at small values of (f> 
and approach a constant value towards (j) = k/ 2. A typical variation of the /-integral 
and «/, distributions for a surface crack loaded under tension is shown in Figure 5-3. 
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Figure 5-3.a Typical J and a /, Distribution along the Surface Crack Front (Tension) 

(a = 4.27 mm, c = 6.97 mm) 

V.D.2 Bending Specimen ./-integral and ai, Distributions 

The stress fields in a surface crack subjected to bending promote a large plastic zone at 
a distance below the free surface of the crack and a much smaller plastic zone towards the 
deepest point of penetration. The stress gradients along the crack front are more severe 
than in tension specimens as evidenced in the characteristic ./-integral and a/, values along 
the crack front. J values reach a maximum value where the plastic zone is largest and 
rapidly decrease with the plastic zone size; however, towards the deepest point of 
penetration constraint calculations for surface crack bend specimens tend to show much 
smaller variations along the crack front and reach a maximum value where the plastic 
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zone is the smallest. A typical surface crack J and a/, plot for a model loaded under 
bending is provided in Figure 5-3. b. 



Figure 5-3.b Typical / and a* Distribution along the Surface Crack Front (Bending) 

(a = 1.516 mm, c = 1.822 mm) 


V.E Fracture Prediction 

The proposed fracture criterion is intended to recognize the location of maximum 
crack extension as the point of highest Ja /, along the crack front. The (f) location of 
maximum Ja/„ denoted as (/) cr n predicted was plotted against the critical 0 location taken 
from the fracture surface images, <j) C nt measured for both the tension and bending cases. 
Figures 5-4. a and 5-4.b show the tension and bending comparisons, respectively. A one- 
to-one correspondence indicates perfect agreement between the measured and predicted 
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critical location. Some test specimens did not exhibit enough crack extension to obtain 
the critical location. These were not included in the comparison. 
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Figure 5-4.a Comparison of Predicted and Measured Critical Location (Bending) 
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Figure 5-4.b Comparison of Predicted and Measured Critical Location (Tension) 


V.F Crack Extension 

A correlation between the crack extension normal to and along the crack front and 
the Joe/, distributions was observed for both the tension and the bending specimens. This 
relationship provides further evidence to the validity of the proposed fracture criterion. 

To observe the relation, the crack extension was digitized and recorded for corresponding 
values of (j) along the crack front. The data was then plotted alongside the Ja/, 
distribution from the appropriate specimen. Sample tension and bending correlation 
examples are presented in Figures 5-5. a and 5-5. b. The cases considered for crack 
extension correlation to Ja/, are given in Appendix C. 
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Figure 5-5.a Sample Crack Extension Correlated with Ja/ t (Bending) 



0.022 

0.020 

0.018 

0.016 

0.014 

0.012 

0.010 

0.008 

0.006 

0.004 

0.002 


Figure 5-5.b Sample Crack Extension Correlated with Ja/ t (Tension) 
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CHAPTER VI 


CONCLUSIONS AND FUTURE WORK 

VI. A Fracture Prediction Validity 

The use of Jan as a fracture prediction criterion for ductile materials is promising. 
For the case of a surface crack under bending, the criterion predicted the critical location 
within approximately 10% error for all specimens considered. For surface cracks under 
tension, the criterion was inconclusive. The results presented for tension showed little 
variation of Joe/, along the crack front in the critical regions, thus a single location of 
maximum Joe/, could not be identified. Thus, crack initiations could occur at 20 / 71 from 
0.25 to 1.0. However, for the tension specimens analyzed, the measured critical location 
fell within the range of constant Joe/, for all considered cases implying that the criterion 
has not been disproved by these results. 

VLB Suggested Future Work 

A rigorous validation of the fracture criterion should be conducted before 
widespread use is considered. Replicate surface cracked specimens of identical crack 
size and applied loading would prove useful in confirming the data presented herein. In 
addition to the application of the criterion to surface cracked geometries, the use of the 
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hyper-local constraint factor as a normalization relating different crack configurations 
would prove invaluable for future development. 
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43 



program ALPHAH 


implicit none 

real::xcoord,ycoord,zcoord,fdestatus,min,ra,rra,distmin,dist,rO,yO,xO,y,x,& 

xdif,ydif,xl,x2,yl,y2,sx,sy,sz,sxy,syz,sxz,Uo,svm,temp,flos,c,a,w,t,l,da_x,& 
tol,aa,bb,cc,counter2,wtsum,ctrad,rad_dist,toler,m 
integer: : i,cmax,bcmax,null 1 ,n,nodeno,NodeMax,j ,k,ir,kk, counter, Elemno.n 1 ,n2,n3 ,n4,& 
n5,n6,n7,n8,ElemMax,ngpt,step,elem,gpt,ii, status 


real,dimension( : , : ),allocatable : :crack, era ckinit,bcrack,bcrackinit,d, crack 1 ,crack2,& 
Node, stress, stravg, alpha 

real,dimension(:),allocatable::GptStrz,GptStrvm,zstress,wtstrz,wt 
integer, dimension! :),allocatable::crackmask,bcmask,pzmask,gptmask 
integer, dimension!:, :),allocatable::cfnmask, element, Const 

!Open Input Files 

open (10, file- crack, crd') 
open (11, fde-belowcrack.crd') 
open (12, file-nodes, crd') 
open (13, file-elements. elm') 

open (14, fde = 'stress.out', STATUS='OLD', ACTION-READ', IOSTAT=status) 

!Open Output Files 

open (20, fde = 'test.out') 
open (21, fde = 'pzone.crd') 
open (22, fde = 'constraint.txt') 

! Input Variables 

flos=1329.7 

c=4.445 lerack length 

a=-3.2 ! crack depth 

w=25.4 ! model width 

t=6.350 Imodel thickness 

1=50.8 Imodel length 

da_x=1.7066660000E-02 Idistance between adjacent nodes in crack tube 
tol=.0009 Itolerance required to select crack front ellipse 

ctrad=0.01 ! crack tube radius 

Unitialize Nodal Coordinate Array to set NodeMax 


n=0 

do while (fdestatus.ge.O) 

read(12,40,iostat=fdestatus) Nodeno,xcoord,ycoord,zcoord 
40 format (18, 3F 17.8) 
n=n+l 
enddo 
n=n-l 

NodeMax=0 



NodeMax=n 

filestatus=0 

rewind(12) 

IBuild Nodal Array 

allocate(Node(NodeMax,4)) 
do i=l,NodeMax 

read(12,*,iostat=fdestatus) Nodeno,xcoord,ycoord,zcoord 
Node(i,l)=Nodeno 
Node(i,2)=xcoord 
Node(i,3)=ycoord 
Node(i,4)=zcoord 
enddo 
filestatus=0 

Unitialize Element Array 
n=0 

do while (filestatns.ge.O) 

read(13,50,iostat=filestatus) Elemno,nl,n2,n3,n4,n5,n6,n7,n8 
50 format (918) 
n=n+l 
enddo 
n=n-l 

ElemMax=0 
ElemMax=n 
filestatus=0 
rewind) 13) 

IBuild Element Array 

allocate(Element(n,9)) 
do i=l,ElemMax 

read) 1 3,50,iostat=fdestatus) Elemno,n 1 ,n2,n3,n4,n5,n6,n7,n8 
Element(i, 1 )=Elemno 
Element(i,2)=nl 
Element(i,3)=n2 
Element(i,4)=n3 
Element(i,5)=n4 
Element(i,6)=n5 
Element(i,7)=n6 
Element(i,8)=n7 
Element(i,9)=n8 
enddo 

Unitialize Crack Array to set Cmax (number of crack front nodes) 


n=0 

do while(fdestatus.ge.O) 
n=n+l 

read) 1 0, *,iostat=filestatus) NULL 1 ,Nodeno,xcoord,ycoord,zcoord 
enddo 
cmax=n- 1 
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f]lestatus=0 

rewind(lO) 

IBuild Initial Crack Array 

allocate(crackinit(cmax,4)) 

30 format(A4,I8,lx,E15.6,lx,E15.6,lx,E15.6) 
do i=l,cmax 

read( 1 0,30,iostat=filestatus) NULL 1 ,Nodeno,xcoord,ycoord,zcoord 
crackinit(i, 1 )=nodeno 
crackinit(i,2)=xcoord 
crackinit(i,3 )=ycoord 
crackinit(i,4)=zcoord 
enddo 
filestatus=0 


Unitialize Below Crack Array to set bcmax (number of nodes on and below the crack front) 


n=0 

do while(filestatus.ge.O) 
n=n+l 

read( 1 1 , *,iostat=filestatus) NULL 1 ,nodeno,xcoord,ycoord,zcoord 
enddo 
filestatus=0 
bcmax=n- 1 
rewind(ll) 

IBuild Initial Below Crack Array 

allocate(bcrackinit(bcmax,4)) 
do i=l, bcmax 

read(ll,30,iostat=filestatus) NULLl,nodeno,xcoord,ycoord,zcoord 
bcrackinit(i,l)=nodeno 
bcrackinit(i,2)=xcoord 
bcrackinit( i, 3 )=ycoord 
bcrackinit(i,4)=zcoord 
enddo 

IFill Crack Front Mask Array and Crack Node/Coordinate Array 

allocate(crackmask(NodeMax)) 
allocate(crack(Nodemax,4)) 
do i=l,cmax 
do j=l,NodeMax 
if (int(crackinit(i,l)).eq.j) then 
crackmask(j)=l 
crack(j,l)=j 

crack(j,2)=crackinit(i,2) 
crack(j ,3 )=crackinit( i, 3 ) 
crack(j,4)=crackinit(i,4) 
endif 
enddo 



enddo 


!Fill Below Crack Mask Array and Node/Coordinate Array 

allocate(bcmask(NodeMax)) 
allocate(bcrack(NodeMax,4)) 
do i=l,bcmax 
do j=l,NodeMax 
if (int(bcrackinit(i,l)).eq.j) then 
bcmask(j)=l 
bcrack(j,l)=j 

bcrack(j,2)=bcrackinit(i,2) 
bcrack(j , 3 )=bcrackinit( i,3 ) 
bcrack(j,4)=bcrackinit(i,4) 
endif 
enddo 
enddo 


lUsing the crack front array which contains all the nodes on the actual crack front 
land the below crack array, which contains the nodes below the crack front, 
Idetermine which node is the shortest distance from each crack front node, 
land store in the array crack 2. The result will be an array of nodes which form a 
Iconcentric ellipse to the crack front 


k=l 

allocate(d(NodeMax,2)) 

allocate(crack2(cmax,3)) 

do i=l,NodeMax 

distmin=99999 

if (crackmask(i).eq.l) then I if node is on crack front 
do j=l,NodeMax 

if (bcmask(j).eq.l) then I if node is below crack front 
if (int(crack(i,l)).ne.int(bcrack(j,l))) then 

dist=sqrt(((bcrack(J,2)-crack(i,2))**2)+((bcrack(j,3)- 

crack(i,3))**2)) 

if (dist.lt.distmin) then 
distmin=dist 
NodeNo=j 
endif 
endif 
endif 
enddo 

d(k,l)=NodeNo 

d(k,2)=distmin 

k=k+l 

endif 

enddo 

k=k-l 


do i=l,cmax 



crack2(i, 1 )=d(i, 1) 
crack2(i,2)=bcrack(d(i,l),2) 
crack2(i,3)=bcrack(d(i,l),3) 
enddo 

allocate(crack 1 (cmax,3 )) 
kk=l 


do i=l,NodeMax 
if (crackmask(i).eq. 1) then 

crack 1 (kk, 1 )=crack(i, 1 ) 
crack 1 (kk,2 )=crack( i,2 ) 
crack 1 (kk,3 )=crack( i, 3 ) 
kk=kk+l 
endif 
enddo 
kk=kk-l 


!call sort_pick(crackl) 

do j=3,k 
aa=crackl(j,3) 
bb=crackl(j,l) 
cc=crackl(j,2) 
do H-1,1,-1 

if (crackl(i,3) <= aa) exit 
crack 1 (i+ 1 , 1 )=crack 1 (i, 1 ) 
crack 1 (i+ 1 ,3 )=crack 1 (i, 3 ) 
crack 1 (i+ 1 ,2)=crack 1 (i,2) 
enddo 

crackl(i+l,3)=aa 

crackl(i+l,l)=bb 

crackl(i+l,2)=cc 

enddo 

!call sort_pick(crack2) 

do j=3,k 
aa=crack2(j,3) 
bb=crack2(j,l) 
cc=crack2(j,2) 
do i=j-l,l,-l 
if (crack2(i,3) <= a) exit 
crack2 (i+ 1 , 1 )=crack2 (i, 1 ) 
crack2 (i+ 1 ,3 )=crack2 (i, 3 ) 
crack2 (i+ 1 ,2)=crack2 (i,2) 
enddo 

crack2(i+l,3)=aa 

crack2(i+l,l)=bb 

crack2(i+l,2)=cc 

enddo 



! Determine which nodes lie in rays 'normal' to the crack front 


allocate(cfnmask(NodeMax,cmax)) 
do i=l,k 

do j=l,NodeMax 
if (bcmask(j).eq. 1) then 
x2=crack2(i,2) 
y2=crack2(i,3) 
xl=crackl(i,2) 
yl=crackl(i,3) 
if (x2-xl.ne.O) then 

X=bcrack(j,2) 

Y=bcrack(j,3) 

m =((y2-y 1 )/ (x2-x 1 )) 
toler=m*(X-x 1 )+y 1 -Y 
toler=abs(toler) 
if (toler.le.da_x) then 
cfnmask(j,i)=l 
else 

cfnmask(j,i)=0 

endif 

else 

Y=bcrack(j,3) 

x=x2-(((y2-Y)*(x2-xl))/(y2-yl)) 
xdif=abs(bcrack(j,2)-x) 
if (xdif.le. 0.00001) then 
cfnmask(j,i)=l 
else 

cfnmask(j,i)=0 

endif 

endif 

endif 

enddo 

enddo 

!Read in results to determine which nodes are yielded 
ngpt=ElemMax* 8 

100 format (Ix,i5,lx,i5,lx,i5,2x,8el4.6) 
do 

read (14,*, IOSTAT=status) temp 
if (status /=0) exit 
enddo 
rewind(14) 

allocate(stress(ngpt,4),STAT=status) 

allocate(gptmask(ngpt)) 

do i= 1 .ngpt 

read(14,100) step,elem,gpt,sx,sy,sz,sxy,syz,sxz,Uo,svm 

stress(i,l)=elem 

stress(i,2)=gpt 



stress(i,3)=sz 

stress(i,4)=svm 

gptmask(i)=l 

enddo 


allocate(GptStrz(NodeMax)) 

allocate(GptStrvm(NodeMax)) 

allocate(StrAvg(NodeMax,3)) 

do j=l,NodeMax 
counter2=0 
do i=l,ngpt 

if (gptmask(i).eq. 1) then 
if (int(stress(i,2)).eq.j) then 
counter2=counter2+ 1 

GptStrz(counter2)=stress(i,3) 

GptStrvm(counter2)=stress(i,4) 

gptmask(i)=0 

endif 

endif 

enddo 

StrAvg(j,l)=j 

StrAvg(j,2)=sum(GptStrz)/counter2 
S tr Avg(j , 3 )=sum( GptS trvm)/ counter2 
do k=l,NodeMax 
GptStrz(k)=0 
GptStrvm(k)=0 
enddo 
enddo 

allocate(Const(NodeMax,3 )) 

do i=l,NodeMax 
Const(i,l)=Node(i,l) 
if (Node(i,4).eq.O) then 

Const(i,3)=l 

else 

Const(i,3)=0 

endif 

enddo 

allocate(pzmask(NodeMax)) 

do i=l,NodeMax 
if (bcmask(i).eq. 1) then 
if (StrAvg(i,3).ge.flos) then 
pzmask(i)=l 
else 

pzmask(i)=0 

endif 

endif 

enddo 



do i=l,NodeMax 
if (pzmask(i).eq. 1) then 
120 format (I8,4F17.8) 

write (21,120) int(Node(i,l)), Node(i,2), Node(i,3), Node(i,4), StrAvg(i,2) 

endif 

enddo 

allocate(zstress(NodeMax)) 

allocate(alpha(cmax,2)) 

allocate(wt(NodeMax)) 

allocate(wtstrz(NodeMax)) 

do i=l,cmax 
counter2=0 
do j=l,NodeMax 
zstress(j)=0 

if (cfnmask(j,i).eq.l) then 
if (pzmask(J).eq. 1) then 
counter2=counter2+ 1 
zstress(j)=StrAvg(j,2) 
endif 
endif 
enddo 

alpha(i,l)=asin(abs(crackl(i,3))/abs(a)) 

alpha(i,2)=(abs(sum(zstress)/counter2))/flos 

enddo 


130 format (3E14.6) 
do i=l,cmax 

write (22,130) alpha(i,l), alpha(i,2) 
enddo 


end ALPHAH 
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Figure B.1-1AT-01 / and ai, Distribution (a = 1.22 mm, c = 1.6 mm, tension) 



Figure B.1-2AT-04 / and ^Distribution (a = 3.20 mm, c = 4.44 mm, tension) 
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Figure B.1-3AT-05 / and «/, Distribution (a = 3.33 mm, c = 4.75 mm, tensioi 



Figure B.1-4AT-06 / and a* Distribution (a = 3.27 mm, c = 4.56 mm, tension) 
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Figure B.1-7AT-09 / and a/, Distribution (a = 4.16 mm, c = 6.63 mm, tension) 



Figure B.1-8BT-01 J and ah Distribution (a = 4.61 mm, c = 9.80 mm, tension) 
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Figure B.1-9CT-01 / and «/, Distribution (a = 3.49 mm, c = 6.44 mm, tension) 



Figure B.l-10 CT-02 / and ^Distribution (a = 3.51 mm, c = 6.65 mm, tension) 
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Figure B.l-11 CT-03 / and a/, Distribution (a = 3.38 mm, c = 6.56 mm, tension) 
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Figure B.l-12 DT-02 / and a* Distribution (a = 1.73 mm, c = 6.42 mm, tension) 
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Figure B.2-1AB-01 / and a i, Distribution (a = 4.356 mm, c = 7.076 mm, bendi 
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Figure B.2-1AB-04 / and a h Distribution (a = 3.185 mm, c = 4.543 mm, bending) 
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Figure B.2-4BB-01 / and a* Distribution (a = 4.62 mm, c = 9.86 mm, bending) 
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Figure B.2-5CB-01 / and a* Distribution (a = 1.712 mm, c = 2.69 mm, bending) 
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Figure B.2-6DB-03 / and ai, Distribution (a = 2.314 mm, c = 6.487 mm, bending) 
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Figure C-l AB-01 Crack Extension Correlated with Ja/, 
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Figure C-2 BB-01 Crack Extension Correlated with Jai, 
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Figure C-3 AT-04 Crack Extension Correlated with Jau 
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Figure C-4 CT-03 Crack Extension Correlated with Jai, 


(} UJ0U o) / m k> p (} UJOU o) / 4 ra p 




